!-------------------------------------------------------------------------------------
subroutine write_derived()
!-------------------------------------------------------------------------------------
  use dimensions
  use parameters

  implicit none

  integer :: fout

  complex    :: epsilon, epsilon0, Csi0, dCsi0dv, dEpsilonDomega
 

  fout = 2

  open(fout, file = trim(runname)//'.out', position='append')

  write(fout,*) 'Derived Parameters:'
  write(fout,*)


  write(fout,*)'   filtered as:  y = (xvec(ix)-Lx/4) / (Lx/03)';
  write(fout,*)'   dkx/kz=', dkx/kz
  write(fout,*)'   actual speckle width = SpeckleWidth x 0.37233 =', &
                     0.37233*SpeckleWidth
  write(fout,*)'   transit time damping rate = u/(actual speckle with) =', &
                     u / (0.37233*SpeckleWidth);
      
         
  ! unit of time is electron plasma frequency of length is electron Debye length
  ! therefore unit of speed on electron thermal speed.
  write(fout,*)'   external potential trapping width = sqrt(phi0) =', sqrt(phi0)
  write(fout,*)'   Domega=', &
                    -(4.366 - 20.86*kz + 28.09* kz**2) * ( exp(-0.5 / kz**2) / kz**2 ), &
                    ' x  phi^0.5 + ', &
                    - 343 * kz**4 * exp( -0.5 / kz**2 ), &
                    ' x phi^1.5'

  write(fout,*)
           
  write(fout,*)'   parallel Courant number =  vPhase * dt / dz =', vPhase*dt/dz
  write(fout,*)'   dt figure of merit for advection = dt kmax Vmax =', &
                     dt*( (nkz/3)-1 )*kz*max(abs(vvec(1)),abs(vvec(jv)))     

  write(fout,*)
  
  write(fout,*)'   max |k| that survives de-Aliasing =',  kxVec(nkx/2  - nkx/6 -1)
  write(fout,"('    number of Fourier modes of a given sign = -1+nkz/3 = ', (I3) )" ) -1+nkz/3

  write(fout,*)

  write(fout,*)"   Dvvvv x pvec(max)^4 x dt =  ", Dvvvv * pvec(jv/2)**4 * dt
  write(fout,*)"   Dvvvv x pvec(max)^4 =       ", Dvvvv * pvec(jv/2)**4
  write(fout,*)"   Dxxxx x kZvec(max)^4 x dt = ", Dxxxx * kZvec(nkz/2)**4 * dt
  write(fout,*)"   Dxxxx x kZvec(max)^4 =      ", Dxxxx * kZvec(nkz/2)**4
  write(fout,*)'   Dxxxx * k_zmax^4 / (k_zmax * dv) =', &
                   Dxxxx * (kZvec(nkz/2  - nkz/6 -1))**3 / dv
  write(fout,*)'                      should be > 1 to resolve 1/kv  resonance'

  write(fout,*)

  write(fout,*)"   kz(max)*v(max) =",  kZvec(nkz/2)*max(abs(vvec(1)),abs(vvec(jv)))
  write(fout,*)"   p(max) = ", pvec(jv/2)
  write(fout,*)


  !-- write "sci" ----------------------------

  call getCsi0(Csi0,cmplx(vPhase,0.0),dCsi0dv)


  write(fout,"('    Re(Csi0) =        ', (e19.8) )" ) real(Csi0)
  write(fout,"('    Im(Csi0/kz^2) =   ', (e19.8) )" ) aimag(Csi0)/kz**2
  write(fout,"('    dCsi0dv =         ', (e19.8), ' + i*'  (e15.8) )" ) &
       real(dCsi0dv),  aimag(dCsi0dv)


  if(  real(Csi0)  >  0  ) then
     epsilon0 = 1 - Csi0/ real(Csi0)
     dEpsilonDomega = - dCsi0dv / sqrt(real(Csi0))**3
     write(fout,"('    kz at resonance = ', (e19.8) )" ) sqrt(real(Csi0))
     write(fout,"('    Landau damping =  ', (e19.8) )" ) &
                  aimag( epsilon0 ) / real( dEpsilonDomega  )
  endif


  epsilon0 = 1 - Csi0/ kz**2

  write(fout,"('    |PhiInternal| =   ', (e19.8) )" ), phi0*abs(   1/epsilon0 - 1 )
  write(fout,*)

  close(fout)


end subroutine write_derived
    

!-------------------------------------------------------------------------------------
  subroutine getCsi0(Csi0,vc,dCsi0dv)     !	real(vc)>=0,  aimag(vc)>=0, or both <0
!-------------------------------------------------------------------------------------

  implicit none
  complex Csi0,vc,dCsi0dv         !	dCsi/dzeta
  real x,mu
  complex Z,Zprime,Zprime2
  
  x=abs ( real(vc) )
  mu=abs( aimag(vc) )
  
  call  getZstuff(x/1.41421,Z,Zprime,Zprime2)     !	1.41421=sqrt(2)
  
  dCsi0dv=3.53553E-01*Zprime2     !	3.53553E-01=1/2**1.5
  
  Csi0=0.5*Zprime + (0.0,1.0)*mu*dCsi0dv
  
  return
  end 


!-------------------------------------------------------------------------------------
  subroutine getZstuff(x,Z,Zprime,Zprime2)        ! should be  good  for  x pos or neg
!-------------------------------------------------------------------------------------
  implicit none
  real x,pi, DAWS, MyDAWS !  DAWS is Dawson integral, as give by IMSL library.
  complex Z,i,Zprime,Zprime2
  data pi,i/3.14159265,(0.0,1.0)/
  
!  Z=i*sqrt(pi)*exp(-x*x) -2*DAWS(x)
   Z=i*sqrt(pi)*exp(-x*x) -2*MyDaws(x)  !  use this if DAWS is unavailable
  
  Zprime=-2*(1+x*Z) !		dZ/dz
  Zprime2=-2*x*Zprime-2*Z     !	d2Z/dx2
  return
  end
  
  real function MyDaws(z1)  !  reproduces DAWS to high accuracy
  implicit none
  real x,z1,Z
  x=abs(z1)
   if( (0.0.le.x).and.(x.lt.1.0) )Z=-0.6953*x**4 + 1.692*x**3  -0.07885*x**2  -1.994*x
   if( (1.0.le.x).and.(x.lt.2.0) )Z=+0.2535577*x**4  -1.905759*x**3 + 5.158011*x**2  -5.463602*x+ 0.8816458
   if( (2.0.le.x).and.(x.lt.2.5) )Z=-(1/x)*( 1+0.5/x**2 +0.75/x**4 +(15./8.)/x**6 ) &
&       +1.563421E-1*x**4  -1.516732*x**3 + 5.520618*x**2  -8.926533*x + 5.400882
   if( x.ge.2.5 )Z=-(1/x)*( 1+0.5/x**2 +0.75/x**4 +(15./8.)/x**6 )  -6.124581*Exp(-3.024965*x)

   MyDaws=-0.5*Z
   if(  z1 < 0  )MyDaws=-MyDaws
   return
   end 
  

!-------------------------------------------------------------------------------------

